{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculation of detectable Gamma ray burst (GRB) event rate associated with Binary Neutron Star (BNS) mergers.\n",
    "\n",
    "* This is without considering a particular detectors' sensitivity, but rather a general calculation of the rate of detectable GRBs based of viewing angle and luminosity distance. I have considered the inclination angle gravitational waves same as the viewing angle wrt GRB jet axis.\n",
    "\n",
    "* For BNS source redshift and mass distribution, refer [here](https://ler.readthedocs.io/en/latest/GW_events.html).\n",
    "\n",
    "* For lens and image properties, refer [here](https://ler.readthedocs.io/en/latest/Lensed_events.html).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equations (GRB)\n",
    "\n",
    "* Ref: Granot, J. & Kumar, P. (2003). \"Distribution of gamma-ray burst ejecta energy with Lorentz factor.\" The Astrophysical Journal, 591(2), 1086-1096. DOI: 10.1086/375575\n",
    "\n",
    "* Disclamer: ChatGPT was used, partially, to get the following expression.\n",
    "\n",
    "* The observed energy of a gamma-ray burst (GRB) as a function of both the angular distance from the jet axis, $\\theta$, and the luminosity distance, $D_L$, is given by:\n",
    "\n",
    "$$E_{\\text{obs}}(\\theta) = \\frac{\\epsilon_0}{4\\pi D_L^2} \\left(\\frac{1}{1 + \\left(\\frac{\\theta}{\\theta_c}\\right)^k}\\right)$$\n",
    "\n",
    "* Here, $E_{\\text{obs}}(\\theta)$ is the observed energy, $\\epsilon_0$ is the energy per unit solid angle at the jet axis, $D_L$ is the luminosity distance, $\\theta$ is the veiwing angle wrt the jet axis, $\\theta_c$ is the core angle of the jet, and $k$ is the power-law index that determines how quickly the energy decreases with angle.\n",
    "\n",
    "* I cosidered the core angle of the jet: $\\theta_c = 5$ deg. (this might be too much)\n",
    "\n",
    "* So, probability of detection is given by:\n",
    "\n",
    "$$P_{det}(\\theta,D_L) = \\left\\{ \n",
    "  \\begin{array}{ c l }\n",
    "    1 & \\text{if } \\theta \\le 5 \\text{ deg and } D_L \\le 46652 \\text{ Mpc } (z\\sim 5) \\\\\n",
    "    1 & \\text{if } \\frac{\\text{a}}{4\\pi D_L^2} \\left(\\frac{1}{1 + \\left(\\frac{\\theta}{5}\\right)^{\\text{c}}}\\right) \\ge 1 \\\\\n",
    "    0 & \\text{otherwise}\n",
    "  \\end{array}\n",
    "\\right.$$\n",
    "\n",
    "* Coefficients (a and c) are solved by considering conditions 2 and 3 below.,\n",
    "  1. If angle <= 5 deg, distance <= distance_cut, then $P_{det}$(bool) = 1.\n",
    "  2. Luminosity(core)/Luminosity(30deg) >= 100; from GW170817.\n",
    "  3. pdet <1, If angle < 30 deg (at distance=40) or distance < 40 Mpc (at angle<30). This is by considering the GRB from GW170817 is barely detectable, with viewing angle=30 deg and distance=40 Mpc. \n",
    "\n",
    "\n",
    "* I have also tried to consider the angular dependence with dipole radiation formula. But, the enenrgy doesn't drop off quick enough in the off-axis viewing angle. So, I have considered the above equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ler.rates import LeR\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Probability of detection of GRBs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finding the coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "c = 2.9542496722537264\n",
      "a = 2513.274122871834\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import fsolve\n",
    "\n",
    "# to find coefficient c\n",
    "# consider: Luminosity(core)/Luminosity(30deg) >= 100; from GW170817\n",
    "def equation(c):\n",
    "    return (1 / (1 + (5 / 5)**c))/(1 / (1 + (30 / 5)**c)) - 100\n",
    "\n",
    "c_guess = 3\n",
    "\n",
    "# Solve the equation\n",
    "c_solution = fsolve(equation, c_guess)\n",
    "\n",
    "print(f\"c = {c_solution[0]}\")\n",
    "\n",
    "# to find coefficient a\n",
    "# consider: pdet = 1, If angle = 30 deg (at distance=40). This is by considering the GRB from GW170817 is barely detectable, with viewing angle=30 deg and distance=40 Mpc.\n",
    "def equation(a):\n",
    "    c = 2.9542496722537264\n",
    "    angle = 30\n",
    "    distance = 40\n",
    "    return (1 / (1 + (angle / 5)**c))* (40/distance)**2 * a/(4*np.pi) - 1\n",
    "\n",
    "a_guess = 3\n",
    "\n",
    "# Solve the equation\n",
    "a_solution = fsolve(equation, a_guess)\n",
    "\n",
    "print(f\"a = {a_solution[0]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function to calculate the probability of detection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "d_L = 46652.22297277721 Mpc\n"
     ]
    }
   ],
   "source": [
    "# Find distance_cut\n",
    "# convert redshift to luminoisty distance\n",
    "import astropy.units as u\n",
    "from astropy.cosmology import LambdaCDM\n",
    "cosmo = LambdaCDM(H0=70, Om0=0.3, Ode0=0.7)\n",
    "z = 5\n",
    "d_L = cosmo.luminosity_distance(z)\n",
    "print(f\"d_L = {d_L}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "# distance_cut: z ~ 5\n",
    "def pdet_grb(angle, distance, distance_cut=46652, duty_cycle=0.5, mask_duty_cycle=True, bool=True):\n",
    "    \"\"\"\n",
    "    Function to calculate the probability of detection for a given angle and distance for GRB. Coefficients are based on 2 and 3 the following conditions,\n",
    "\n",
    "    0. GRB jet, core angle <= 5 deg\n",
    "    1. If angle <= 5 deg, distance <= distance_cut, then pdet(bool) = 1\n",
    "    2. Luminosity(core)/Luminosity(30deg) >= 100; from GW170817\n",
    "    3. pdet <1, If angle < 30 deg (at distance=40) or distance < 40 Mpc (at angle<30). This is by considering the GRB from GW170817 is barely detectable, with viewing angle=30 deg and distance=40 Mpc.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    angle : numpy.ndarray\n",
    "        Angle between the GRB jet and viewing angle in (rad)\n",
    "    distance : numpy.ndarray\n",
    "        Distance between the GRB source and the satellite in (Mpc)\n",
    "    distance_cut : float\n",
    "        Core angular size of the GRB jet in (rad)\n",
    "        default is 25422.742 Mpc\n",
    "    duty_cycle : float\n",
    "        Duty cycle of detector(s)\n",
    "        default is 0.5\n",
    "    bool : bool\n",
    "        if True, return absolute value of pdet\n",
    "        if False, return the pdet value as boolean (with duty cycle applied)\n",
    "    \"\"\"\n",
    "\n",
    "    # coefficients, obtained with scipy.optimize fsolve by considering the condition 2,3 listed in docstring\n",
    "    c = 2.9542496722537264\n",
    "    a = 2513.274122871834\n",
    "\n",
    "    # convert angle to degree\n",
    "    angle = np.degrees(angle)\n",
    "\n",
    "    # make sure that the input data is a numpy array\n",
    "    angle, distance = np.array([angle]).reshape(-1), np.array([distance]).reshape(-1)\n",
    "\n",
    "    # angle should be less than 90 or equal to 90\n",
    "    if len(angle[angle > 90]) > 0:\n",
    "        angle[angle > 90] = angle[angle > 90]%90\n",
    "\n",
    "    if bool:\n",
    "        # calculate the probability of detection, absolute value\n",
    "        pdet = abs((1 / (1 + (angle / 5)**c))* (40/distance)**2 * a/(4*np.pi))\n",
    "\n",
    "        # find idx of angle <= 5 and distance <= distance_cut\n",
    "        idx = (angle <= 5) & (distance <= distance_cut)\n",
    "        # apply the condition, condition 1 from docstring\n",
    "        pdet[idx] = 1\n",
    "\n",
    "        if mask_duty_cycle:\n",
    "            # apply the duty cycle\n",
    "            # sample random numbers from 0 to 1 and check if it is less than the duty cycle\n",
    "            num_ = np.random.rand(len(angle))\n",
    "            mask_duty_cycle = num_ > duty_cycle\n",
    "            pdet[mask_duty_cycle] = 0\n",
    "\n",
    "        # return the pdet as boolean\n",
    "        return (pdet>=1).astype(int)\n",
    "    else:\n",
    "        # return the probability of detection (absolute value)\n",
    "        return abs((1 / (1 + (angle / 5)**c))* (40/distance)**2 * a/(4*np.pi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 0 1 0]\n"
     ]
    }
   ],
   "source": [
    "angle = np.array([1, 1, 30, 90])\n",
    "angle = np.radians(angle)\n",
    "distance = np.array([46652, 46653, 40, 40])\n",
    "print(pdet_grb(angle, distance, mask_duty_cycle=False, bool=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pdet condition checks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "condition 1, pdet(angle=core_angle,distance=distance_cut):  [1]\n",
      "condition 2, absolute value pdet(angle=5,distance=40)/pdet(angle=30,distance=40):  [100.]\n",
      "condition 3\n",
      "  i) Detectable, pdet(angle=30,distance=40):  [1]\n",
      "  ii) Not-Detectable, pdet(angle=31,distance=40):  [0]\n",
      "  iii) Not-Detectable, pdet(angle=30,distance=50):  [0]\n"
     ]
    }
   ],
   "source": [
    "print(\"condition 1, pdet(angle=core_angle,distance=distance_cut): \", pdet_grb(angle=np.radians(5), distance=46652, mask_duty_cycle=False, bool=True))\n",
    "\n",
    "print(\"condition 2, absolute value pdet(angle=5,distance=40)/pdet(angle=30,distance=40): \", pdet_grb(angle=np.radians(5), distance=40, mask_duty_cycle=False, bool=False)/pdet_grb(angle=np.radians(30), distance=40, mask_duty_cycle=False, bool=False))\n",
    "\n",
    "print(\"condition 3\")\n",
    "print(\"  i) Detectable, pdet(angle=30,distance=40): \", pdet_grb(angle=np.radians(30), distance=40, mask_duty_cycle=False, bool=True))\n",
    "print(\"  ii) Not-Detectable, pdet(angle=31,distance=40): \", pdet_grb(angle=np.radians(31), distance=40, mask_duty_cycle=False, bool=True))\n",
    "print(\"  iii) Not-Detectable, pdet(angle=30,distance=50): \", pdet_grb(angle=np.radians(30), distance=50, mask_duty_cycle=False, bool=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* let's write out the function so that it can be used in with LeR\n",
    "* consider 50% duty cycle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gwsnr.utils import save_json_dict\n",
    "\n",
    "# let's write out the function so that it can be used in with LeR\n",
    "def pdet_calculator(gw_param_dict, duty_cycle=0.5, mask_duty_cycle=True, bool=True, output_jsonfile=False):\n",
    "    \"\"\"\n",
    "    Function to calculate the probability of detection for a given angle and distance for GRB. This is based on the following condition\n",
    "\n",
    "    1. GRB jet, core angle <= 5 deg\n",
    "    2. If angle <= 5 deg, distance <= distance_cut, then pdet = 1\n",
    "    3. Luminosity(core)/Luminosity(30deg) > 100, from GW170817 \n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    gw_param_dict : dict\n",
    "        dictionary containing the parameters for the GW event\n",
    "    \"\"\"\n",
    "\n",
    "    # get the angle and distance from the dictionary\n",
    "    angle = gw_param_dict['theta_jn']\n",
    "    distance = gw_param_dict['luminosity_distance']\n",
    "\n",
    "    # calculate the probability of detection\n",
    "    pdet = pdet_grb(angle, distance, duty_cycle=duty_cycle, mask_duty_cycle=mask_duty_cycle, bool=bool)\n",
    "    pdet_net_dict = dict(pdet_net=pdet)\n",
    "\n",
    "    # Save as JSON file, if output_jsonfile is provided\n",
    "    if output_jsonfile:\n",
    "        output_filename = (\n",
    "            output_jsonfile if isinstance(output_jsonfile, str) else \"pdet.json\"\n",
    "        )\n",
    "        save_json_dict(pdet_net_dict, output_filename)\n",
    "\n",
    "    # return the pdet\n",
    "    return pdet_net_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'pdet_net': array([1, 1, 0])}\n"
     ]
    }
   ],
   "source": [
    "# test\n",
    "gw_param_dict = {'theta_jn': np.radians(np.array([1, 30, 90])), 'luminosity_distance': np.array([2500, 40, 40])}\n",
    "print(pdet_calculator(gw_param_dict, mask_duty_cycle=False, bool=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rate calculation of GRB with LeR\n",
    "\n",
    "* initialize LeR with pdet calculator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "ler = LeR(verbose=False, pdet_finder=pdet_calculator, event_type='BNS')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Unlensed events"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "unlensed params will be store in ./ler_data/unlensed_param_grb.json\n",
      "chosen batch size = 100000 with total size = 500000\n",
      "There will be 5 batche(s)\n",
      "Batch no. 1\n",
      "sampling gw source params...\n",
      "calculating pdet...\n",
      "Batch no. 2\n",
      "sampling gw source params...\n",
      "calculating pdet...\n",
      "Batch no. 3\n",
      "sampling gw source params...\n",
      "calculating pdet...\n",
      "Batch no. 4\n",
      "sampling gw source params...\n",
      "calculating pdet...\n",
      "Batch no. 5\n",
      "sampling gw source params...\n",
      "calculating pdet...\n",
      "saving all unlensed_params in ./ler_data/unlensed_param_grb.json...\n"
     ]
    }
   ],
   "source": [
    "ler.batch_size = 100000\n",
    "unlensed_param = ler.unlensed_cbc_statistics(size=500000, output_jsonfile='unlensed_param_grb.json', resume=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lensed params will be store in ./ler_data/lensed_param_grb.json\n",
      "chosen batch size = 100000 with total size = 500000\n",
      "There will be 5 batche(s)\n",
      "Batch no. 1\n",
      "sampling lensed params...\n",
      "solving lens equations...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████| 100000/100000 [00:25<00:00, 3996.84it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Invalid sample found. Resampling 4 lensed events...\n",
      "solving lens equations...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████| 4/4 [00:02<00:00,  1.38it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "calculating pdet...\n",
      "Batch no. 2\n",
      "sampling lensed params...\n",
      "solving lens equations...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████| 100000/100000 [00:24<00:00, 4073.71it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Invalid sample found. Resampling 3 lensed events...\n",
      "solving lens equations...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████| 3/3 [00:02<00:00,  1.09it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "calculating pdet...\n",
      "Batch no. 3\n",
      "sampling lensed params...\n",
      "solving lens equations...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████| 100000/100000 [00:24<00:00, 4064.70it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Invalid sample found. Resampling 2 lensed events...\n",
      "solving lens equations...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████| 2/2 [00:02<00:00,  1.35s/it]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "calculating pdet...\n",
      "Batch no. 4\n",
      "sampling lensed params...\n",
      "solving lens equations...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████| 100000/100000 [00:24<00:00, 4016.73it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "calculating pdet...\n",
      "Batch no. 5\n",
      "sampling lensed params...\n",
      "solving lens equations...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████| 100000/100000 [00:25<00:00, 3942.88it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Invalid sample found. Resampling 2 lensed events...\n",
      "solving lens equations...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████| 2/2 [00:02<00:00,  1.35s/it]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "calculating pdet...\n",
      "saving all lensed_params in ./ler_data/lensed_param_grb.json...\n"
     ]
    }
   ],
   "source": [
    "ler.batch_size = 100000\n",
    "lensed_param = ler.lensed_cbc_statistics(size=500000, output_jsonfile='lensed_param_grb.json', resume=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "getting unlensed_params from json file ./ler_data/unlensed_param_grb.json...\n",
      "given detectability_condition == 'pdet'\n",
      "total unlensed rate (yr^-1) (with step function): 2109.796632137801\n",
      "number of simulated unlensed detectable events: 10190\n",
      "number of all simulated unlensed events: 500000\n",
      "storing detectable unlensed params in ./ler_data/unlensed_param_detectable.json\n",
      "getting lensed_params from json file ./ler_data/lensed_param_grb.json...\n",
      "given detectability_condition == 'pdet'\n",
      "total lensed rate (yr^-1) (with pdet function): 2.117598781203888\n",
      "number of simulated lensed detectable events: 4261\n",
      "number of simulated all lensed events: 500000\n",
      "storing detectable lensed params in ./ler_data/lensed_param_detectable.json\n",
      "unlensed_rate (per year): 2109.796632137801\n",
      "lensed_rate (per year): 2.117598781203888\n",
      "ratio: 996.3155678330854\n"
     ]
    }
   ],
   "source": [
    "rate_ratio, unlensed_param_detectable, lensed_param_detectable = ler.rate_comparison_with_rate_calculation(\n",
    "    unlensed_param='unlensed_param_grb.json',\n",
    "    lensed_param='lensed_param_grb.json',\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* we have 2110 GRBs (detectable, un-lensed) per year\n",
    "* we have 2.11 GRBs (detectable, lensed) per year\n",
    "* ratio of lensed to un-lensed is 1:1000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rate calculation of GW with LeR\n",
    "\n",
    "* using internal snr_calculator function, with [gwsnr](https://gwsnr.readthedocs.io/en/latest/) package.\n",
    "* Sensitivity: O4 design sensitivity, [L1,H1,V1]\n",
    "* for detection criteria of un-lensed BNS GW refer [here](https://ler.readthedocs.io/en/latest/GW_events.html), and for lensed BNS GW refer [here](https://ler.readthedocs.io/en/latest/Lensed_events.html). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "ler = LeR(verbose=False, event_type='BNS')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "unlensed params will be store in ./ler_data/unlensed_param_bns_gw.json\n",
      "chosen batch size = 100000 with total size = 500000\n",
      "There will be 5 batche(s)\n",
      "Batch no. 1\n",
      "sampling gw source params...\n",
      "calculating snrs...\n",
      "Batch no. 2\n",
      "sampling gw source params...\n",
      "calculating snrs...\n",
      "Batch no. 3\n",
      "sampling gw source params...\n",
      "calculating snrs...\n",
      "Batch no. 4\n",
      "sampling gw source params...\n",
      "calculating snrs...\n",
      "Batch no. 5\n",
      "sampling gw source params...\n",
      "calculating snrs...\n",
      "saving all unlensed_params in ./ler_data/unlensed_param_bns_gw.json...\n"
     ]
    }
   ],
   "source": [
    "ler.batch_size = 100000\n",
    "unlensed_param = ler.unlensed_cbc_statistics(size=500000, output_jsonfile='unlensed_param_bns_gw.json', resume=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lensed params will be store in ./ler_data/lensed_param_bns_gw.json\n",
      "chosen batch size = 100000 with total size = 500000\n",
      "There will be 5 batche(s)\n",
      "Batch no. 1\n",
      "sampling lensed params...\n",
      "solving lens equations...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████| 100000/100000 [00:25<00:00, 3861.70it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "calculating snrs...\n",
      "Batch no. 2\n",
      "sampling lensed params...\n",
      "solving lens equations...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████| 100000/100000 [00:24<00:00, 4054.86it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Invalid sample found. Resampling 2 lensed events...\n",
      "solving lens equations...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████| 2/2 [00:02<00:00,  1.37s/it]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "calculating snrs...\n",
      "Batch no. 3\n",
      "sampling lensed params...\n",
      "solving lens equations...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████| 100000/100000 [00:24<00:00, 4052.37it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Invalid sample found. Resampling 2 lensed events...\n",
      "solving lens equations...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████| 2/2 [00:02<00:00,  1.45s/it]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "calculating snrs...\n",
      "Batch no. 4\n",
      "sampling lensed params...\n",
      "solving lens equations...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████| 100000/100000 [00:24<00:00, 4021.48it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Invalid sample found. Resampling 1 lensed events...\n",
      "solving lens equations...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████| 1/1 [00:03<00:00,  3.19s/it]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "calculating snrs...\n",
      "Batch no. 5\n",
      "sampling lensed params...\n",
      "solving lens equations...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████| 100000/100000 [00:24<00:00, 4062.46it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Invalid sample found. Resampling 2 lensed events...\n",
      "solving lens equations...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████| 2/2 [00:03<00:00,  1.52s/it]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "calculating snrs...\n",
      "saving all lensed_params in ./ler_data/lensed_param_bns_gw.json...\n"
     ]
    }
   ],
   "source": [
    "ler.batch_size = 100000\n",
    "lensed_param = ler.lensed_cbc_statistics(size=500000, output_jsonfile='lensed_param_bns_gw.json', resume=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "getting unlensed_params from json file ./ler_data/unlensed_param_bns_gw.json...\n",
      "given detectability_condition == 'step_function'\n",
      "total unlensed rate (yr^-1) (with step function): 2.898641104016606\n",
      "number of simulated unlensed detectable events: 14\n",
      "number of all simulated unlensed events: 500000\n",
      "storing detectable unlensed params in ./ler_data/unlensed_param_detectable.json\n",
      "getting lensed_params from json file ./ler_data/lensed_param_bns_gw.json...\n",
      "given detectability_condition == 'step_function'\n",
      "total lensed rate (yr^-1) (with step function): 0.008448528345568198\n",
      "number of simulated lensed detectable events: 17\n",
      "number of simulated all lensed events: 500000\n",
      "storing detectable lensed params in ./ler_data/lensed_param_detectable.json\n",
      "unlensed_rate (per year): 2.898641104016606\n",
      "lensed_rate (per year): 0.008448528345568198\n",
      "ratio: 343.094203564711\n"
     ]
    }
   ],
   "source": [
    "rate_ratio, unlensed_param_detectable, lensed_param_detectable = ler.rate_comparison_with_rate_calculation(\n",
    "    unlensed_param='unlensed_param_bns_gw.json',\n",
    "    lensed_param='lensed_param_bns_gw.json',\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* we have 2.9 GWs (detectable, un-lensed) per year\n",
    "* we have 0.0084 GWs (detectable, lensed) per year\n",
    "* ratio of lensed to un-lensed is 1:343"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Final results\n",
    "\n",
    "### Table for GRB and GW rates\n",
    "\n",
    "| Event | Unlensed | Lensed | ratio |\n",
    "| --- | --- | --- | --- |\n",
    "| GRB | 2110 | 2.11 | 1:1000 |\n",
    "| GW | 2.9 | 0.0084 | 1:343 |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ler",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}